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Abstract 

The B-spline basis set method is applied to determining the rovibrational eigen-spectrum of 
diatomic molecules. A particular attention is paid to a challenging numerical task of an accurate 
and efficient description of the vibrational levels near the dissociation limit (halo-state and Feshbach 
molecules). Advantages of using B-splines are highlighted by comparing the performance of the 
method with that of the commonly-used discrete variable representation (DVR) approach. Several 
model cases, including the Morse potential and realistic potentials with 1/R? and l/i? 6 long-range 
dependence of the internuclear separation are studied. We find that the B-spline method is superior 
to the DVR approach and it is robust enough to properly describe the Feshbach molecules. The 
developed numerical method is applied to studying the universal relation of the energy of the last 
bound state to the scattering length. We numerically illustrate the validity of the quantum-defect- 
theoretic formulation of such a relation for a 1/-R 6 potential. 
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I. INTRODUCTION 



Finite basis sets technique is an important numerical tool in solving quantum mechanical 
problems, e.g., in quantum chemistry One of the popular recent developments is the 
use of B-splines in such calculations. In atomic physics, the applications of B-splines were 
stimulated by Walter R. Johnson's work [2] and here, in this special issue dedicated to 
celebrating his contributions to atomic physics, we are delighted to present yet another 
robust application of B-splines. 

The reason to the popularity of the B-splines in practical applications is due to the 
fact that they form a sufficiently complete basis set with a reasonably small number of 
basis functions. Numerical accuracy of the calculations approaches that of the traditional 
finite-difference methods, such as the Numerov method 3], with the advantage of a global 
non-iterative determination of eigenergies and eigenstates. 

Here we apply the B-spline method to study rovibrational eigen-spectrum of diatomic 
molecules and compare the performance of the method with that of the discrete variable 
representation (DVR) approach. Previously, the B-spline method was successfully applied 
to finding vibrational spectrum of the Morse potential in Refs. {3, 4|. Here we focus on 
the more challenging problem of describing vibrational states near the dissociation limit of 
realistic long-range potentials. One difficulty lies in the variation of the local de Broglie 
wavelength by several orders of magnitude from the short range to the long range region. 
Several authors jf], have discussed how the efficiency of the DVR methods could be 
improved via the implementation of a mapping procedure, where the grid step is adapted to 
the variation of the de Broglie wavelength. In the present paper we compare the mapped sine 
grid method of Ref. Q to B-splines calculations also using a mapping procedure. Molecular 
bound states near dissociation limit play an important role in the formation of ultracold 
molecules ?], Q] and in the determination of scattering properties in the low-energy regime, 
in particular scattering length or more generally the threshold energy-dependence of the 
phaseshifts. The vibrational wavefunctions then extend to distances much larger than the 
typical length of the chemical bond. Recently several experimental groups succeeded in 
making loosely-bound ultracold molecules by sweeping B-fields through the magnetically- 



induced Feshbach resonances (see review 



10J). Such Feshbach molecules may be considered 



as halo-state systems, since the vibrational wavefunctions extend well into the classically-for- 
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bidden region. While the halo-state systems, due to their universal behavior for a wide range 



of quantum- mechanical systems, deserve a special attention on their own right [llj, there 
are emerging applications based on the Feshbach molecules: for example, several schemes 

n n 

of transferring Feshbach molecules to lower vibrational levels and down to v — [12|, [13J 
have been proposed. In the prerequisite numerical time-dependent studies, an expansion 
over a suitably-chosen quasi-spectrum is required, and the initial state near dissociation 
limit has to be well represented by this quasi-spectrum. The challenge there is the accurate 
representation of the evanescent part of the wavefunction in the non classical region. B- 
splines, with their superior numerical performance, demonstrated here, may prove useful 
in such theoretical studies. We shall therefore evaluate this performance by comparing to 
analytical results when available (bound levels of the Morse potential) or to well-established 
numerical methods. 

Motivated by the spectacular developments in low-energy collision physics of ultracold 
atoms, the universal laws governing near-threshold physics have generated a considerable 
interest over the last decade. In particular, here, with the developed numerical method, 
we investigate a relation of the energy of the last bound state to the scattering length. 
For potentials without a long range tail, such a relation is a well-known prediction of the 
effective-range theory (see e.g., 14|). For the van der Waals potentials, the effective range 
theory has to be improved to account for their asymptotic behavior and Gao [15] has recently 
derived the proper law relating these two quantities in the framework of the quantum defect 
theory (QDT). Here, using the developed B-spline code, we verify numerically the validity 
of this new formulation. We find that compared to the effective-range result, the QDT 
expression remains accurate over a much wider range of parameters than expected. 

The paper is organized as follows. First, in Section [Til we set-up the numerical method 
using the Galerkin technique and expansion of the molecular wavefunctions over the B-spline 
basis. We also describe an efficient molecular grid used in the calculations and recapitu- 
late main features of the DVR method. In Section IHIl we apply the method to finding 
ro-vibrational spectra of various potentials and compare the results with those from the 
DVR method. We start with the Morse potential, where analytical results are available, 
and proceed to realistic potentials, varying with the internuclear separations, R, as 1/R 3 
and 1/R 6 at large R. With the developed method, in Section ITV] we analyze the relation 
between the scattering length and the position of the last bound state and compare our 
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numerical results with the predictions of the QDT and the effective-range theories. Finally, 
the conclusions are drawn in Section [V] 

II. PROBLEM SETUP 

We are interested in solving the radial time-independent Schrodinger equation for vibra- 
tional motion of nuclei of a diatomic molecule 

__L„; M + ^ M + Zg±l))„ j(Ji) . 

Euj(R), (1) 

where [i is the reduced molecular mass, J is the rotational quantum number and V (R) 
is the electronic Born-Oppenheimer potential. Unless specified otherwise, atomic units, 
h = \e\ = m e = 1, are used throughout. 

A. B-spline approach 

General mathematical introduction to B-splines and a collection of codes to manipulate 
these basis functions may be found in Ref . [JiJ . Here we briefly recapitulate properties of the 
B-splines relevant to our discussion. We deal with a set of n functions defined on a support 
grid {ti}. A B-spline, B^ (R) , number i of order k is a piecewise polynomial of degree k — 1 
inside an interval of the support grid ti < R < t i+ k. It vanishes outside this support interval. 
The B-splines are positive functions on their support interval. In applications, the common 
choice (also used here) is to make the end-points of the support grid fc-fold degenerate, 

t\ = ti — ■ ■ ■ — tk — Rmim 

where n is the total number of B-splines in the set. With such a choice of the grid, the 
first B-spline, (R) is the only spline which does not vanish at -R m i n . Similarly, the only 
non- vanishing B-spline at the end-point R max is the last B-spline, B^ n (R). Notice that the 
spline support grid ti directly maps on the radial grid, except the multiply defined end-points 
that map onto the first and last points of the radial grid. 
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Below, we employ the Galerkin method to obtain a quasi-spectrum of the radial 



Schrodinger equation, see, e.g., 



171 ] . Central to this approach is an observation that the 



differential equation ([T]) may be derived by seeking an extremum of the action integral, 



R„ 

-E / u 2 j(R)dR. 



nun 



Further, we expand the ro-vibrational wavefunctions in terms of the B-spline set, 

n-l 

uj (R) = £ dB? } (R) ■ (2) 

t=2 

Notice that we discard the first and the last B-spline of the set to enforce the boundary 
conditions uj (-R m i n ) = and uj (-R max ) = 0. The remaining splines vanish identically at the 
end-points of the grid. 

We substitute the expansion ([2]) in the action integral and seek its extremum with respect 
to the expansion coefficients. As a result, we arrive at the generalized eigen- value equation 
for the vector of the coefficients c = (c 2 , c 3 , ...e n _i): 



with matrices 



Ac = E B c , (3) 



/•Rmax r 1 dB . dB . 

13 ~ J Rmin UlRtlR + 

2B i (v{R) + ^±^)BA ( in. (1) 



2R 2 



/"-Rmax 

B{j = / B{Bj dR. 



-^min 

The resulting eigenf unctions are orthonormal and form a numerically complete basis set in 
the space of piece-wise polynomials of order k—1. The choice of the number of basis functions 
is determined by the nodal structure of the wavefunctions that we wish to represent. 



B. Mapped grid method 

Choosing numerical grid for solving the radial Schrodinger equation for loosely bound 
molecules requires special consideration. Realistic potentials support a large number of 
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'nstead, we employ a more 
5), [f| , first implemented in 



bound states. Near the dissociation limit the corresponding wavefunctions have a large 
number of nodes. Moreover, the distance between two nodes, and hence the local De Broglie 
wavelength, grows larger as we approach the outer turning point of the potential. A large 
fraction of the wavefunction (especially for halo-state molecules) may reside in the classically- 
forbidden region. Because of this behavior of the vibrational states, here we depart from the 
usual choice of the radial grid of a constant step as in Refs. 3,0]. 
efficient grid as prescribed by the "mapped grid" method of Ref. 
the framework of the DVR method described below. 

In the "mapped grid" method the radial grid is based on the adaptive coordinate defined 

as 

x (R) = /T 1 ^ f dR'^E max -V cnv (R'), 

where V env (R) is the enveloping potential (it is chosen to be either the same as or slightly 
deeper than the original potential V(R)), R m i n is somewhat smaller than the position of 
the repulsive inner part of the potential, E max is the maximum energy for which accurate 
results are wanted and p max is the corresponding value of the total linear momentum. The 
grid transformation x (R) efficiently rescales the radial coordinate by the local de Broglie 
wavelength. Factor (3 < 1 makes the radial step smaller than the local de Broglie wavelength 
and improves the representation of the wavefunction in the classically-forbidden region. We 
use a constant step of Ax = vr^,/p max for the adaptive coordinate. This choice translates 
into a variable step of the radial grid, 

At this point we recast the solution of the differential equation in terms of the generalized 
eigenvalue equation (j3J). To solve this problem, we developed a numerical code using B-spline 



routines of Ref. 



• 3- 



Below we evaluate the performance of the method by studying the 



rovibrational spectrum of several potentials. 



C. Discrete variable approach 



The DVR approach to the computation of vibrational wavefunctions [18j, is based on 
a collocation scheme. A wavefunction ip is approximated by its projection Pip on a linear 
combination of iV interpolation functions, such that <p and Pip have the same values at the 



collocation points. The wavefunction ip is thus represented by its values at the collocation 
points. The Hamiltonian is represented by a matrix, which can be used to compute bound 
and continuum states or to simulate the temporal evolution of a wavepacket. Spectral and 



collocation methods are discussed in a famous monograph by D. Gottlieb and S. Orszag [19]. 

A great variety of systems have been studied, using various sets of orthogonal interpola- 
tion functions. In contrast with the B-splines, such functions do not vanish outside a small 
interval, but rather they all are defined on the whole grid, and differ by the number of nodes. 

For applications to ultracold molecules, with bound and quasi-bound vibrational levels 
in asymptotically R~ 6 and R~ 3 potentials, Kokoouline et al l2l| have implemented 

a Mapped Fourier Grid method where the interpolation functions are plane waves. The 
grid step is rescaled to the value of the local de Broglie wavelength, as described above 
in Eq.Q. Accurate results were obtained both for the vibrational energies and for the 
wavef unctions, using a number of basis functions slightly larger than the number of nodes 
of the wavefunction of the upper level . The accuracy could be checked by comparison with 
asymptotic methods 22] derived from generalized quantum defect theory . However, the 
occurrence of ghosts levels after diagonalization of the Hamiltonian matrix appeared as a 
drawback of the mapping procedure. Willner et al [6| have shown that when replacing the 
plane waves by a basis of N sine functions of the adaptive coordinate x, 

8 k (x) = J^em(kj-x > ) (k = l,...,N-l), (6) 

with nodes at both ends of the grid, most of the ghost levels would disappear. 
The relevant formulae for the collocation scheme can be found in Ref. 6j. Note that the 
number of basis functions is entirely determined by the number of grid steps. The length of 
the grid is related to the constant grid step 5x in the x coordinate by 

L = N5x (7) 

Levels of the Cs2 dimer with a binding energy as small as ~ 10~ 16 a. u. could be computed, 
for which the vibrational wavefunction extends up to 100 000 ao, i.e. a few tens of microns. 
This wavefunction with 528 nodes is computed with a grid of only 706 points: it is typical of 
a halo molecule, most of the probability density lying in the classically forbidden region. The 
efficiency of a set of oscillating sine functions to represent this slowly decreasing exponential 
function is then questionable. A discussion on the appearance of ghost levels shows that 



they are influenced by the value chosen for the parameter (3 : a compromise has to be 
found between the suppression of ghost states (/3 small) and a minimum value of grid points 
{(3 ~ 1). Moreover, the numerical representation of the potential, where an analytical long 
range behavior is usually matched to an interpolation function between discrete ab initio 
data at short range may be a source of unphysical levels Our choice in the present paper is 
to compare the efficiency of the B-spline and sine-grid methods for the same grid. 



III. NUMERICAL EXAMPLES 



A. Morse potential 



As a test of the quality of our numerical approach, we start with the Morse potential |23| , 
which has no long range tail but has an advantage of having analytically known energy levels 
and wavef unctions. The Morse potential is given by 

V (r) = D ( e - 2a(r - ro) - 2e- a(r " ro) ) , (8) 

where D is the dissociation energy, ro is the equilibrium position, and the parameter a 
governs the spatial extent of the potential. The energies of the bound states are known 
exactly, 



E v = —D + fko (v + 1/2) - \ hw (v + l/2) z , (9) 

where the vibrational quantum number v = 0,1, ..Vd, with the maximum, vd = 
[a~ 1 ^2fxD — 1/2J . In these formulas, the vibrational frequency is 

[2D\ 1/2 , . 

= • do) 

In calculations we use Morse potential fitted to the ground state potential of 133 CS2 
dimer. The parameters of the employed Morse potential are (in atomic units) r = 8.77, 
D = 0.016627, a = 0.372031199. This potential supports 170 bound states. 

We carry out the DVR and B-spline computations using identical grids. Given the same 
grid, the accuracy of the resulting eigen-values depends only on the basis, sin (DVR) or 
B-spline set, and the method of solution of the Schrodinger equation (collocation versus 
Galerkin method). In Table [J, we compare the computed energies (both DVR and B- 
splines) with analytical results for vibrational levels near the dissociation limit. The results 
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marked a were computed using a relatively small grid of N = 275 points (i? m in = 6.3 do, 
-Rmax = 100 a , and ft = 0.7). The larger and denser grid (entries marked 6) has iV = 553 
points, -R m i n = 6.3 do, -Rmax = 2000 d , and ft = 0.4. In both cases E max = 10~ 8 . The order 
of B- splines is k — 15. 

First we consider a case of the coarse grid (a). The accuracy of reproducing the energies 
of the low-lying states in the B-spline method is at the level of 10~ n cm -1 , while the DVR 
method has an accuracy of about 10~ 7 cm -1 . More substantial is the difference in the 
spectrum near the dissociation limit. Here the DVR spectrum is perturbed by a "ghost" 
state v = 168. Because of the ghost state, the resulting number of bound states in the DVR 
spectrum is incorrect. The spectral position of the "ghost" state varies as the parameters 
of the grid change; for example, the bound spectrum is no longer perturbed in case of the 
larger grid (b). By contrast, the B-spline set spectrum is free of the ghost states regardless 
of the choice of the grid. 

As we shift to the denser grids (case (b)), the numerical accuracy of both methods im- 
proves. Because of the improved accuracy, in Table [I] we list deviations of the numerical 
energies from the analytical values. Again, we observe that the B-spline method outper- 
forms the DVR method in terms of accuracy. This conclusion seem to hold irrespective 
of a particular choice of grid. The accuracy of computing the energy of the last bound 
level requires special consideration. The relevant outer classical turning point is located at 
R = 50.6 do- However, the wavefunction substantially extends into the classically forbidden 
region. The small grid (-R m ax = 100 do) can not fully accommodate this tail. As the size 
of the cavity is increased to 2000 do for the large grid (b), the B-spline method starts to 
recover 4-5 significant figures of the exact result for the energy of the last bound state. Yet 
the DVR method reproduces only the leading significant figure. 

The superior performance of the B-spline method seems to be due to the compactness of 
B-splines. A given B-spline extends only over k intervals of the grid: the B-spline number 
i vanishes identically outside a support interval (ij,tj + &). In particular, it means that for 
a given coordinate R only a sum of k (in our case k = 15) B-splines contributes. This is 
in a stark contrast to the DVR method: here all N ~ 1000 rapidly oscillating functions 
contribute to a value of the wavefunction at a given coordinate, leading to the deterioration 
of numerical accuracy. Moreover, it is intuitively clear that while the DVR sin basis is natural 
for describing rapid oscillations in the classically-allowed region, the forbidden region with 
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V 


Analytical 


E v , DVR a 


E v , B-splines a 


AE V , DVR 6 


AE V B-splines b 


162 


-8.2264075 


-8.2263504 


-8.2264050 


q v in- 


-6 


5 x 10" 13 


163 


-6.3205792 


-6.3205024 


-6.3205766 


i v in- 

± s\ _L\J 


-5 


-1 x 10- 12 


164 


-4.6655178 


-4.6654299 


-4.6655114 


i v in- 


-5 


1 x 10~ 12 


165 


-3.2612233 


-3.2611065 


-3.2612136 


o v in- 

L A _L\J 


-5 


8 x 10" 11 


166 


-2.1076957 


-2.1075487 


-2.1074619 


2 x 10" 


-5 


9 x 10" 9 


167 


-1.2049349 


-1.2046878 


-1.2023402 


4 x 10" 


-5 


2 x 10" 6 


168 


-0.5529410 


-1.0604859 


-0.54813941 


6 x 10" 


-5 


9 x 10~ 6 


169 


-0.15171398 


-0.5527958 


-0.14926508 


2 x 10" 


-4 


4 x 10" 6 


170 


-1.2538365 x 10~ 3 


-0.1507383 


-1.0319455 x 10" 


" 3 2 x 10- 


-4 


9 x 10~ 8 



Table I: Comparison of the accuracy of the DVR and B-spline methods in the case of the Morse 
potential for two choices of radial grids. Results marked (a) are for the case of a coarse grid and 
results marked (b) are for a finer grid. 



its extended exponential tail requires well-balanced interference of many basis functions. 
The accurate description of the classically-forbidden region becomes more important as we 
approach the dissociation limit. Namely in this limit the advantages of using B-splines 
become more substantial. 



B. Attractive l/R 3 interactions 

Compared to the Morse potential, realistic molecular potential display a long range 
tail leading to a dense vibrational spectrum near the dissociation limit. The long-range 
neutral-atom interactions depend on the internuclear distance as —C n /R n , with n > 3. 

The most challenging is the case of two atoms interacting via attractive —C3/R 3 , C3 > 0, 
interactions. Such potentials, for example, do not possess scattering length [24]. As a 
particular example, we consider the A l Y^ potential of 87 Rb2 dimer correlating to the 5s + 5p 
asymptotic limit, shown in the upper panel of Fig. [TJ This potential is attractive at large 
internuclear distances, V (R) ~ —C 3 /R 3 . In our specific case C 3 ~ 17.81 a.u.. 



As shown by Le Roy and Bernstein 22| for long-range potentials varying as V (R) 
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-C 3 /R* 

E V = D-[H 3 (v D -v)}\ (11) 

where the constant if 3 is related to the long-range constant. In our case the dissociation 
limit D = 0. 

We plot our computed dependence of (— E v ) l / & on the vibrational quantum number in 
the lower panel of Fig. [TJ We see that the Le Roy-Bernstein formula, Eq. pip , is followed 
up to v ~ 435. This equation was derived u sing semi-classical arguments and it is known 
to be violated for the last vibrational levels 25]. However, in our case the deviation from 
Eq. fllip for levels of v > 435 is simply due to limitations of the double precision arithmetics 
(15 significant figures) used in the computations. Indeed, the energy spectrum spans 14 
orders of magnitude: the lowest vibrational state has an energy of —2.9 x 10~ 2 a.u., while 
E v= 435 ~ —3.8 x 1CT 16 . Both the B-spline and the DVR methods, since they reproduce the 
entire spectrum in one shot, do not cope well with the loss of numerical accuracy. If desired, 
numerical accuracy could be improved by switching to quadruple precision arithmetics. 

We find that the B-spline results for levels v < 435 were insensitive to a particular choice 
of the grid, as long as the -R max wa s well beyond the outer classical turning point of the 
wavefunction. By contrast, the DVR code has produced a multitude of ghost levels, and, 
for the best choice of the grid parameters, we were able to reproduce positions of at most 
430 vibrational levels. 

The computed wavefunction of the v = 434 level is plotted in Fig. [2l For this state, the 
classical turning point is located at 1.9 x 10 5 bohr. The B-spline code was run using the 
mapped grid with the following parameters: -R m i n = 5.0 a , -R ma x = 1 x 10 7 a , (3 = 0.5, 
-Emax = 10~ 15 . This corresponds to 1292 grid-points. Notice that the v = 434 wavefunction 
has 434 nodes, yet it was accurately computed using only 1292 grid-points. This is an 
excellent demonstration of the efficiency of the mapped-grid technique coupled with the 
B-spline method. 



IV. RELATION BETWEEN THE POSITION OF THE LAST BOUND LEVEL 
AND THE SCATTERING LENGTH 

Here we consider two atoms interacting at long-range separations via attractive —Cq/R & , 
Cq > 0, potentials. We will employ two scaling parameters: the van der Waals length 
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Figure 1: Upper panel: Molecular potential A Sj of Rb2 molecule. Lower panel: comparison of 
the Le Roy-Bernstein fit (solid line) with the results obtained with the B-spline code (squares). 

0.004 
0.003 
0.002 

X 

£ 0.001 

a 



a 

# 0.000 
-0.001 
-0.002 

1 y 1 : 2x10 s 3x10 5 4x10" 5x10 5 

R, bohr 

Figure 2: Vibrational wavefunction of v = 434, J = level of the Rb A 1 T,^ electronic potential 
as computed in the B-spline method. The vertical line marks the position of the classical turning 
point. 

f 6 = (2yuC 6 ) 1//4 and energy E$ = 1/ (fi (f 6 ) 2 ). In particular, the regime of quantum halo 
states is reached when the energy of the last bound state is — E_i <^ E 6 and its spatial 
extension reaches distances much larger than f 6 . 
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We have investigated the performance of the B-spline method in the case of a realistic 
molecular potential that follows the 1/R 6 power law at large distances (this is the case of 
the ground state of the alkali dimers). The numerical results are quite similar to the already 
presented cases of the Morse and 1/R? long-range potential. Instead, in this section we 
use the developed method to study universal relation between the scattering length and the 
position of the last bound state in the molecular potential. To this end we focus on a simple 
model of hard-core sphere with the van der Waals tail. In this model the short-range physics 
is modeled by placing an impenetrable wall at R = R : 

, x oo , R < Rq t . 

V(R) ={ . (12) 

( -C 6 /R 6 , R > R 

This simple model offers insights into the universal laws of low-energy scattering. Let us 
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enumerate several analytical results 
are formulated in terms of the scaling factor 

2tt 



26j for this model relevant to our discussion. These 



(r(i/4)) 2 

and accumulated phase inside the potential 



f 6 « 0.477989 f 6 , 



$ - 6 



2*8 



which determines the physics close to threshold. Indeed, the number of bound states is given 



by 



26] 



N b = L$/tt - 7/8J + 1 , 



and the scattering length a by 26| 



a = a (1 - tan($ - 3tt/8)) . (13) 

For | a | /re ^> 1, there is either a bound level close to the dissociation limit (a > 0) or a 
virtual state (a < 0). 



In our numerical study, we take Cq = 6851 for the ground-state Cs dimer [27J], and a 
reduced mass for 133 Cs atoms. For 133 Cs2 molecule f 6 « 202 a , and E e pa 4.4 x 10 -5 cm -1 . 
Increasing Rq, the position of the inner "hard" wall of the potential, reduces the number 
of bound states in the potential. For example, we find from analytical formula that a new 
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bound state appears at the value of -Rq ~ 6.02073 ao- The potential binds 180 states for Rq 
just below Rq and 179 states for R just above Rq. 

For our initial numerical test we choose the position of the inner wall at Rq = 6.02 bohr. 
B-spline method reliably produces all 179 bound states and reveals a loosely bound state 
with the energy of —9.33 x 10~ 12 a.u. We verified that the energies of the states near the 
dissociation limit follow the Le Roy-Bernstein pattern (similar to the analysis presented in 
Fig. [T]for the l/R? potential.) In this case, however, some additional observations can be 
made. 

For Rq = 6.02 ao, the scattering length, Eq. ( fl3|) is large and positive, a = +796 a . Large 
and positive scattering lengths result from having a bound state just below the threshold. 
In this regime, the energy of the last bound state may be approximated by [ijj] 

E?W « _jL !_ x 
_1 2/j(a-a) 2 

l + Cl -^- + c 2r ^V (14) 
[a — a) {a — a) z J 

where c\ ~ 0.4387552, c 2 ~ —0.2163139. The above expression was derived using the 
quantum defect theory and it substantially differs from the commonly-used effective range 
expansion formula 

From Eq. (1141) . we find E^ T « -9.35 x 10~ 12 , while the effective -range formula results 
in E_i ~ —6.52 x 10 -12 . Clearly, our numerical result, —9.33 x 10~ 12 a.u., supports the 



analytical analysis 



151 ] . In this calculation, the parameters of the grid were chosen to be 



-Rmin = Ro, -Rmax = 5 x 10 4 , (3 = 0.4, with the number of points 2407. When the number 
of points was reduced by a factor of 3, the energy of the last bound state was affected in 
the third significant figure. We again notice that the DVR method was unable to match the 
numerical accuracy of the B-spline approach. 

While offering an improved accuracy over the effective-range expression, the QDT Eq. (fT4"j) 
is still an approximate result. In Fig. [3], we compare the QDT prediction with our numerical 
results. Here we move the position of the inner wall just below the critical value of Rq ~ 
6.02073, at which the least bound state disappears. The range of the values for the position 
of the inner wall was chosen so that the scattering length remained positive. An increase 
in Rq translates into increasingly larger values of the scattering length. For a/a ^> 1, i.e., 
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near the threshold, both the effective-range and the QDT results become identical. As the 
scattering length decreases, the effective range approximation rapidly loses its accuracy. Our 
comparison in Fig. [3] clearly demonstrates that, compared to the conventional effective-range 
theory, the QDT approximation is applicable over a much wider range of parameters. At 
the same time, as Ro is decreased from its critical value, the QDT approximation starts to 
break down at Rq ~ 6.017 bohr. The relevant parameter governing the validity of Eq. (1141) is 
the reduced scattering length a/a: the critical value Rq ~ 6.017 ao corresponds to a/a w 2. 
To reiterate, the QDT formula, Eq. fUl) . is an excellent approximation as long as a/ a > 2, 
while the effective range approximation requires a/ a 3> 1. 

Finally, it is worth pointing out that our method is robust enough to reproduce halo 
states of diatomic molecules bound by the van der Waals forces. We varied R just below the 
threshold value and examined the energies of the least bound state produced by the B-spline 
method. For example, for Rq = 6.0207, we obtain with the B-spline code E_i = —1.71 x 
10~ 14 a.u., while analytical results are E^ T « -1.82 x 10" 14 a.u. and w -1.78 x 10~ 14 
a.u. The binding energies are four orders of magnitude smaller than the van der Waals 
energy. At the same time, the corresponding scattering length, governing the extent of the 
wavefunction, is about 2 x 10 4 bohr, i.e., two orders of magnitude larger than the van der 
Waals length. Satisfying both enumerated conditions signifies reaching the universal regime 
of quantum halo states. 



V. CONCLUSION 



With the experimental control of quantum-mechanical systems becoming more refined, 
new theoretical tools have to be adopted to meet the new challenges. Recently, frag- 
ile Feshbach (quantum halo-state) molecules became an experimental reality (see e.g., 



Ref. 



28|). Motivated by this progress, here we developed a numerical method for solv- 



ing the Schrodinger equation for diatomic molecules based on the B-spline finite basis sets. 
The method produces a numerically complete quasi-spectrum of rovibrational states. We 
find, that B-splines offer an accurate description of the loosely-bound molecular states near 
the dissociation limit. The quasi-spectrum is entirely devoid of the unphysical ghost states 



which appear in DVR method and require special effort to be eliminated 



coupled with the "mapped grid" method of Ref. Q], the representation is both accurate and 



Moreover, 
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Figure 3: Energy of the last bound state as a function of the position of the inner wall of the model 
potential. Dots mark numerical results obtained with the B-spline approach. Predictions of the 
effective-range approximation are shown with a dashed line and that of the quantum-defect theory 
- with a solid line. 

efficient: both rapidly-oscillating part of the wavefunction in the classically-allowed region 
and the slowly-varying exponential tail in the classically-forbidden region are adequately 
reproduced. As an application of the developed method we investigated the universal law 
relating the energy of the last bound state to the scattering length. We find that the new 



QDT formulation of such a law for 1/R G potentials by Gao 15J remains valid over a substan- 



tially wider range of parameters than the commonly-used effective-range approximation. 
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